Beta diversity
Paths and libraries setting
# Load main packages, paths and custom functions
source("../../../source/main_packages.R")
source("../../../source/paths.R")
source("../../../source/functions.R")
# Load supplementary packages
packages <- c("openxlsx", "kableExtra", "dynamicTreeCut", "cowplot", "ggpubr")
invisible(lapply(packages, require, character.only = TRUE))HCA on raw data
Preparation
# load rdata file with phyloseq object
setwd(path_rdata)
ps <- readRDS("MED_phyloseq.rds")
# transform count table into percent table
ps.percent <- transform_sample_counts(ps, (function(x) x / sum(x)))
# extract otu and tax tables
otu <- data.frame(ps.percent@otu_table)
tax <- data.frame(ps.percent@tax_table)
sam <- data.frame(ps.percent@sam_data)
# add Genus column in otu table
otu$Genus <- tax$Genus
otu <- otu %>% select(c(Genus, everything()))
# select Wolbachia
otu_wolbachia <- otu[otu$Genus=="Wolbachia" & !is.na(otu$Genus),] %>% t()
otu_wolbachia <- otu_wolbachia[-1,]
# add sample names as rownames
samples_names <- rownames(otu_wolbachia)
# convert column into numeric type
otu_wolbachia <- data.frame(apply(otu_wolbachia, 2, function(x) as.numeric(as.character(x))))
# add column Wolbachia with sum of Wolbachia percentage
otu_wolbachia$Wolbachia <- rowSums(otu_wolbachia)
# set rownames
rownames(otu_wolbachia) <- samples_names
# add column Sample with sample id
otu_wolbachia$Sample <- rownames(otu_wolbachia)
otu_wolbachia <- otu_wolbachia %>% select(c(Sample, Wolbachia))
# print
otu_wolbachia %>%
kbl(booktable=TRUE) %>%
kable_paper("hover", full_width = F)| Sample | Wolbachia | |
|---|---|---|
| CTC1 | CTC1 | 0.0000000 |
| CTC10 | CTC10 | 0.0352626 |
| CTC11 | CTC11 | 0.0077123 |
| CTC12 | CTC12 | 0.0191972 |
| CTC13 | CTC13 | 0.0026615 |
| CTC14 | CTC14 | 0.1074139 |
| CTC15 | CTC15 | 0.0094047 |
| CTC2 | CTC2 | 0.0087316 |
| CTC3 | CTC3 | 0.0007014 |
| CTC4 | CTC4 | 0.0088827 |
| CTC5 | CTC5 | 0.0017734 |
| CTC6 | CTC6 | 0.0020249 |
| CTC7 | CTC7 | 0.0050378 |
| CTC9 | CTC9 | 0.0058233 |
| NP14 | NP14 | 0.9888611 |
| NP2 | NP2 | 0.9993923 |
| NP20 | NP20 | 0.4705882 |
| NP27 | NP27 | 0.9789303 |
| NP29 | NP29 | 0.0492611 |
| NP30 | NP30 | 0.2894737 |
| NP34 | NP34 | 0.3157895 |
| NP35 | NP35 | 0.0032998 |
| NP36 | NP36 | 0.0000000 |
| NP37 | NP37 | 0.0002480 |
| NP38 | NP38 | 0.0218636 |
| NP39 | NP39 | 0.0023176 |
| NP41 | NP41 | 0.0000162 |
| NP42 | NP42 | 0.0004591 |
| NP43 | NP43 | 0.0070779 |
| NP44 | NP44 | 0.0000127 |
| NP5 | NP5 | 0.9999973 |
| NP8 | NP8 | 0.9999492 |
| S100 | S100 | 1.0000000 |
| S101 | S101 | 0.8433759 |
| S102 | S102 | 1.0000000 |
| S104 | S104 | 0.9999237 |
| S105 | S105 | 0.9999100 |
| S106 | S106 | 1.0000000 |
| S107 | S107 | 0.9999808 |
| S108 | S108 | 0.9999641 |
| S109 | S109 | 0.9999492 |
| S110 | S110 | 1.0000000 |
| S113 | S113 | 0.9941796 |
| S121 | S121 | 1.0000000 |
| S122 | S122 | 1.0000000 |
| S123 | S123 | 1.0000000 |
| S124 | S124 | 0.9999449 |
| S126 | S126 | 0.9997374 |
| S127 | S127 | 1.0000000 |
| S128 | S128 | 1.0000000 |
| S146 | S146 | 0.9993204 |
| S147 | S147 | 0.9999205 |
| S148 | S148 | 0.9993646 |
| S150 | S150 | 1.0000000 |
| S151 | S151 | 1.0000000 |
| S152 | S152 | 0.9998674 |
| S153 | S153 | 0.9998826 |
| S154 | S154 | 0.9993767 |
| S160 | S160 | 0.9999035 |
| S162 | S162 | 0.9999603 |
| S163 | S163 | 0.9999188 |
| S164 | S164 | 0.9993741 |
| S165 | S165 | 0.9995488 |
| S166 | S166 | 0.9992847 |
| S167 | S167 | 0.9997153 |
| S169 | S169 | 1.0000000 |
| S170 | S170 | 1.0000000 |
| S18 | S18 | 0.9909091 |
| S19 | S19 | 0.6928881 |
| S20 | S20 | 0.9376633 |
| S21 | S21 | 0.6642608 |
| S22 | S22 | 0.5062214 |
| S23 | S23 | 0.9421566 |
| S24 | S24 | 0.9054080 |
| S25 | S25 | 0.9970647 |
| S26 | S26 | 0.9733037 |
| S27 | S27 | 0.2155884 |
| S28 | S28 | 0.9646902 |
| S29 | S29 | 0.0098804 |
| S30 | S30 | 0.1662131 |
| S31 | S31 | 0.1847797 |
| S32 | S32 | 0.2115416 |
| S33 | S33 | 0.3796025 |
| S34 | S34 | 0.5882353 |
| S35 | S35 | 0.4507416 |
| S36 | S36 | 0.1426408 |
| S37 | S37 | 0.0420938 |
| S38 | S38 | 0.8371053 |
| S39 | S39 | 0.3274431 |
| S40 | S40 | 0.4336303 |
| S41 | S41 | 0.3043478 |
| S42 | S42 | 0.9997565 |
| S43 | S43 | 0.9369888 |
| S44 | S44 | 0.9795664 |
| S45 | S45 | 0.9815447 |
| S46 | S46 | 0.4921260 |
| S47 | S47 | 0.4081812 |
| S48 | S48 | 0.9938841 |
| S49 | S49 | 0.9925667 |
| S50 | S50 | 0.9898880 |
| S51 | S51 | 0.9936881 |
| S52 | S52 | 0.9833898 |
| S53 | S53 | 1.0000000 |
| S55 | S55 | 0.9999604 |
| S56 | S56 | 1.0000000 |
| S57 | S57 | 0.9998373 |
| S58 | S58 | 0.9977978 |
| S59 | S59 | 1.0000000 |
| S60 | S60 | 1.0000000 |
| S61 | S61 | 1.0000000 |
| S63 | S63 | 0.9955108 |
| S64 | S64 | 1.0000000 |
| S65 | S65 | 0.3235294 |
| S77 | S77 | 0.0000000 |
| S79 | S79 | 1.0000000 |
| S80 | S80 | 1.0000000 |
| S83 | S83 | 1.0000000 |
| S84 | S84 | 1.0000000 |
| S85 | S85 | 0.9994486 |
| S86 | S86 | 1.0000000 |
| S87 | S87 | 0.9998787 |
| S88 | S88 | 0.9977787 |
| S89 | S89 | 0.0032180 |
HCA
Run HCA
Dendogram
# 3 groups
plot(hca.ward, ask=FALSE)
rect.hclust(hca.ward, k=3)
# set displaying
rh <- rect.hclust(hca.ward, k=3, border = c("red", "blue", "grey"))beg_clus <- head(cumsum(c(1, lengths(rh))), -1)
beg_clus <- beg_clus + 5
y_clus <- weighted.mean(rev(hca.ward$height)[2:3], c(4, 1))
# plot
plot(hca.ward, cex=0.7)
rect.hclust(hca.ward, k=3, border = c("red", "blue", "grey"))
text(x=beg_clus, y=y_clus, col=c("red","blue", "grey"), labels=c("High infection", "Low infection", "Medium infection"), font=2)# save plot
setwd(path_plot)
tiff("1Ec_dendogram.tiff", units="in", width=18, height=8, res=300)
plot(hca.ward, cex=0.7, xlab="Samples", sub="")
rect.hclust(hca.ward, k=3, border = c("red", "blue", "grey"))
text(x=beg_clus, y=y_clus, col=c("red","blue", "grey"), labels=c("High infection", "Low infection", "Medium infection"), font=2)
dev.off()## quartz_off_screen
## 2
png("1Ec_dendogram.png", units="in", width=18, height=8, res=300)
plot(hca.ward, cex=0.7, xlab="Samples", sub="")
rect.hclust(hca.ward, k=3, border = c("red", "blue", "grey"))
text(x=beg_clus, y=y_clus, col=c("red","blue", "grey"), labels=c("High infection", "Low infection", "Medium infection"), font=2)
dev.off()## quartz_off_screen
## 2
Split
## CTC1 CTC10 CTC11 CTC12 CTC13 CTC15 CTC2 CTC3 CTC4 CTC5 CTC6 CTC7 CTC9
## 1 1 1 1 1 1 1 1 1 1 1 1 1
## NP29 NP35 NP36 NP37 NP38 NP39 NP41 NP42 NP43 NP44 S29 S37 S77
## 1 1 1 1 1 1 1 1 1 1 1 1 1
## S89 CTC14 NP20 NP30 NP34 S19 S21 S22 S27 S30 S31 S32 S33
## 1 2 2 2 2 2 2 2 2 2 2 2 2
## S34 S35 S36 S39 S40 S41 S46 S47 S65 NP14 NP2 NP27 NP5
## 2 2 2 2 2 2 2 2 2 3 3 3 3
## NP8 S100 S101 S102 S104 S105 S106 S107 S108 S109 S110 S113 S121
## 3 3 3 3 3 3 3 3 3 3 3 3 3
## S122 S123 S124 S126 S127 S128 S146 S147 S148 S150 S151 S152 S153
## 3 3 3 3 3 3 3 3 3 3 3 3 3
## S154 S160 S162 S163 S164 S165 S166 S167 S169 S170 S18 S20 S23
## 3 3 3 3 3 3 3 3 3 3 3 3 3
## S24 S25 S26 S28 S38 S42 S43 S44 S45 S48 S49 S50 S51
## 3 3 3 3 3 3 3 3 3 3 3 3 3
## S52 S53 S55 S56 S57 S58 S59 S60 S61 S63 S64 S79 S80
## 3 3 3 3 3 3 3 3 3 3 3 3 3
## S83 S84 S85 S86 S87 S88
## 3 3 3 3 3 3
Make a dataframe from these groups
## CTC1 CTC10 CTC11 CTC12 CTC13 CTC15 CTC2 CTC3 CTC4 CTC5 CTC6 CTC7 CTC9
## 1 1 1 1 1 1 1 1 1 1 1 1 1
## NP29 NP35 NP36 NP37 NP38 NP39 NP41 NP42 NP43 NP44 S29 S37 S77
## 1 1 1 1 1 1 1 1 1 1 1 1 1
## S89 CTC14 NP20 NP30 NP34 S19 S21 S22 S27 S30 S31 S32 S33
## 1 2 2 2 2 2 2 2 2 2 2 2 2
## S34 S35 S36 S39 S40 S41 S46 S47 S65 NP14 NP2 NP27 NP5
## 2 2 2 2 2 2 2 2 2 3 3 3 3
## NP8 S100 S101 S102 S104 S105 S106 S107 S108 S109 S110 S113 S121
## 3 3 3 3 3 3 3 3 3 3 3 3 3
## S122 S123 S124 S126 S127 S128 S146 S147 S148 S150 S151 S152 S153
## 3 3 3 3 3 3 3 3 3 3 3 3 3
## S154 S160 S162 S163 S164 S165 S166 S167 S169 S170 S18 S20 S23
## 3 3 3 3 3 3 3 3 3 3 3 3 3
## S24 S25 S26 S28 S38 S42 S43 S44 S45 S48 S49 S50 S51
## 3 3 3 3 3 3 3 3 3 3 3 3 3
## S52 S53 S55 S56 S57 S58 S59 S60 S61 S63 S64 S79 S80
## 3 3 3 3 3 3 3 3 3 3 3 3 3
## S83 S84 S85 S86 S87 S88
## 3 3 3 3 3 3
# change colnames
colnames(wolbachia.groups) <- "Group"
# change rownames
wolbachia.groups$Sample <- rownames(wolbachia.groups)
# add Group to otu_wolbachia with merge
wolbachia.groups <- wolbachia.groups %>% merge(otu_wolbachia %>% select(c(Sample, Wolbachia)), by="Sample")
# change rownames
rownames(wolbachia.groups) <- wolbachia.groups$Sample
# print
wolbachia.groups %>%
kbl(booktable=TRUE) %>%
kable_paper("hover", full_width = F)| Sample | Group | Wolbachia | |
|---|---|---|---|
| CTC1 | CTC1 | 1 | 0.0000000 |
| CTC10 | CTC10 | 1 | 0.0352626 |
| CTC11 | CTC11 | 1 | 0.0077123 |
| CTC12 | CTC12 | 1 | 0.0191972 |
| CTC13 | CTC13 | 1 | 0.0026615 |
| CTC14 | CTC14 | 2 | 0.1074139 |
| CTC15 | CTC15 | 1 | 0.0094047 |
| CTC2 | CTC2 | 1 | 0.0087316 |
| CTC3 | CTC3 | 1 | 0.0007014 |
| CTC4 | CTC4 | 1 | 0.0088827 |
| CTC5 | CTC5 | 1 | 0.0017734 |
| CTC6 | CTC6 | 1 | 0.0020249 |
| CTC7 | CTC7 | 1 | 0.0050378 |
| CTC9 | CTC9 | 1 | 0.0058233 |
| NP14 | NP14 | 3 | 0.9888611 |
| NP2 | NP2 | 3 | 0.9993923 |
| NP20 | NP20 | 2 | 0.4705882 |
| NP27 | NP27 | 3 | 0.9789303 |
| NP29 | NP29 | 1 | 0.0492611 |
| NP30 | NP30 | 2 | 0.2894737 |
| NP34 | NP34 | 2 | 0.3157895 |
| NP35 | NP35 | 1 | 0.0032998 |
| NP36 | NP36 | 1 | 0.0000000 |
| NP37 | NP37 | 1 | 0.0002480 |
| NP38 | NP38 | 1 | 0.0218636 |
| NP39 | NP39 | 1 | 0.0023176 |
| NP41 | NP41 | 1 | 0.0000162 |
| NP42 | NP42 | 1 | 0.0004591 |
| NP43 | NP43 | 1 | 0.0070779 |
| NP44 | NP44 | 1 | 0.0000127 |
| NP5 | NP5 | 3 | 0.9999973 |
| NP8 | NP8 | 3 | 0.9999492 |
| S100 | S100 | 3 | 1.0000000 |
| S101 | S101 | 3 | 0.8433759 |
| S102 | S102 | 3 | 1.0000000 |
| S104 | S104 | 3 | 0.9999237 |
| S105 | S105 | 3 | 0.9999100 |
| S106 | S106 | 3 | 1.0000000 |
| S107 | S107 | 3 | 0.9999808 |
| S108 | S108 | 3 | 0.9999641 |
| S109 | S109 | 3 | 0.9999492 |
| S110 | S110 | 3 | 1.0000000 |
| S113 | S113 | 3 | 0.9941796 |
| S121 | S121 | 3 | 1.0000000 |
| S122 | S122 | 3 | 1.0000000 |
| S123 | S123 | 3 | 1.0000000 |
| S124 | S124 | 3 | 0.9999449 |
| S126 | S126 | 3 | 0.9997374 |
| S127 | S127 | 3 | 1.0000000 |
| S128 | S128 | 3 | 1.0000000 |
| S146 | S146 | 3 | 0.9993204 |
| S147 | S147 | 3 | 0.9999205 |
| S148 | S148 | 3 | 0.9993646 |
| S150 | S150 | 3 | 1.0000000 |
| S151 | S151 | 3 | 1.0000000 |
| S152 | S152 | 3 | 0.9998674 |
| S153 | S153 | 3 | 0.9998826 |
| S154 | S154 | 3 | 0.9993767 |
| S160 | S160 | 3 | 0.9999035 |
| S162 | S162 | 3 | 0.9999603 |
| S163 | S163 | 3 | 0.9999188 |
| S164 | S164 | 3 | 0.9993741 |
| S165 | S165 | 3 | 0.9995488 |
| S166 | S166 | 3 | 0.9992847 |
| S167 | S167 | 3 | 0.9997153 |
| S169 | S169 | 3 | 1.0000000 |
| S170 | S170 | 3 | 1.0000000 |
| S18 | S18 | 3 | 0.9909091 |
| S19 | S19 | 2 | 0.6928881 |
| S20 | S20 | 3 | 0.9376633 |
| S21 | S21 | 2 | 0.6642608 |
| S22 | S22 | 2 | 0.5062214 |
| S23 | S23 | 3 | 0.9421566 |
| S24 | S24 | 3 | 0.9054080 |
| S25 | S25 | 3 | 0.9970647 |
| S26 | S26 | 3 | 0.9733037 |
| S27 | S27 | 2 | 0.2155884 |
| S28 | S28 | 3 | 0.9646902 |
| S29 | S29 | 1 | 0.0098804 |
| S30 | S30 | 2 | 0.1662131 |
| S31 | S31 | 2 | 0.1847797 |
| S32 | S32 | 2 | 0.2115416 |
| S33 | S33 | 2 | 0.3796025 |
| S34 | S34 | 2 | 0.5882353 |
| S35 | S35 | 2 | 0.4507416 |
| S36 | S36 | 2 | 0.1426408 |
| S37 | S37 | 1 | 0.0420938 |
| S38 | S38 | 3 | 0.8371053 |
| S39 | S39 | 2 | 0.3274431 |
| S40 | S40 | 2 | 0.4336303 |
| S41 | S41 | 2 | 0.3043478 |
| S42 | S42 | 3 | 0.9997565 |
| S43 | S43 | 3 | 0.9369888 |
| S44 | S44 | 3 | 0.9795664 |
| S45 | S45 | 3 | 0.9815447 |
| S46 | S46 | 2 | 0.4921260 |
| S47 | S47 | 2 | 0.4081812 |
| S48 | S48 | 3 | 0.9938841 |
| S49 | S49 | 3 | 0.9925667 |
| S50 | S50 | 3 | 0.9898880 |
| S51 | S51 | 3 | 0.9936881 |
| S52 | S52 | 3 | 0.9833898 |
| S53 | S53 | 3 | 1.0000000 |
| S55 | S55 | 3 | 0.9999604 |
| S56 | S56 | 3 | 1.0000000 |
| S57 | S57 | 3 | 0.9998373 |
| S58 | S58 | 3 | 0.9977978 |
| S59 | S59 | 3 | 1.0000000 |
| S60 | S60 | 3 | 1.0000000 |
| S61 | S61 | 3 | 1.0000000 |
| S63 | S63 | 3 | 0.9955108 |
| S64 | S64 | 3 | 1.0000000 |
| S65 | S65 | 2 | 0.3235294 |
| S77 | S77 | 1 | 0.0000000 |
| S79 | S79 | 3 | 1.0000000 |
| S80 | S80 | 3 | 1.0000000 |
| S83 | S83 | 3 | 1.0000000 |
| S84 | S84 | 3 | 1.0000000 |
| S85 | S85 | 3 | 0.9994486 |
| S86 | S86 | 3 | 1.0000000 |
| S87 | S87 | 3 | 0.9998787 |
| S88 | S88 | 3 | 0.9977787 |
| S89 | S89 | 1 | 0.0032180 |
Groups statistics
With raw labels
# convert summary to table
summary.table <- function(df, n, var){
var <- df[df$Group==n, var]
x <- summary(var)
a <- data.frame(x=matrix(x),row.names=names(x))
colnames(a) <- n
return(a)
}
# make table with stats
wolbachia.summary <- summary.table(wolbachia.groups, 1, "Wolbachia") %>% cbind(summary.table(wolbachia.groups, 2, "Wolbachia")) %>% cbind(summary.table(wolbachia.groups, 3, "Wolbachia"))
# print
wolbachia.summary %>%
kbl(booktable=TRUE) %>%
kable_paper("hover", full_width = F)| 1 | 2 | 3 | |
|---|---|---|---|
| Min. | 0.0000000 | 0.1074139 | 0.8371053 |
| 1st Qu. | 0.0005802 | 0.2155884 | 0.9948452 |
| Median | 0.0032998 | 0.3274431 | 0.9998826 |
| Mean | 0.0091467 | 0.3654875 | 0.9891909 |
| 3rd Qu. | 0.0091437 | 0.4705882 | 1.0000000 |
| Max. | 0.0492611 | 0.6928881 | 1.0000000 |
With new labels
## [1] "1" "2" "3"
# change levels of Group
wolbachia.groups[wolbachia.groups$Group==1, "Group"] <- "Low infection"
wolbachia.groups[wolbachia.groups$Group==2, "Group"] <- "Medium infection"
wolbachia.groups[wolbachia.groups$Group==3, "Group"] <- "High infection"
# make table with stats
wolbachia.summary.labels <- summary.table(wolbachia.groups, "Low infection", "Wolbachia") %>% cbind(summary.table(wolbachia.groups, "Medium infection", "Wolbachia")) %>% cbind(summary.table(wolbachia.groups, "High infection", "Wolbachia"))
# print
wolbachia.summary.labels %>%
kbl(booktable=TRUE) %>%
kable_paper("hover", full_width = F)| Low infection | Medium infection | High infection | |
|---|---|---|---|
| Min. | 0.0000000 | 0.1074139 | 0.8371053 |
| 1st Qu. | 0.0005802 | 0.2155884 | 0.9948452 |
| Median | 0.0032998 | 0.3274431 | 0.9998826 |
| Mean | 0.0091467 | 0.3654875 | 0.9891909 |
| 3rd Qu. | 0.0091437 | 0.4705882 | 1.0000000 |
| Max. | 0.0492611 | 0.6928881 | 1.0000000 |
HCA on qPCR data
Preparation
# load qPCR file
setwd(path_qpcr)
qpcr <- read.xlsx("quanttiWolbachiaschriekeetal.xlsx")
# select Target/ref column
qpcr <- qpcr[,c("Sample.Name", "Target/Ref")]
colnames(qpcr) <- c("Sample", "Ratio_Target_Ref")
# remove the bad S45 sample
qpcr <- qpcr[-28,]
rownames(qpcr) <- qpcr$Sample
# samples without Wolbachia
qpcr_no_wolbachia <- qpcr[qpcr$Ratio_Target_Ref==0,]
qpcr_no_wolbachia$Group <- 0
qpcr_no_wolbachia <- qpcr_no_wolbachia %>% select(-c(Ratio_Target_Ref))
# samples with Wolbachia
qpcr_wolbachia <- qpcr[qpcr$Ratio_Target_Ref!=0,]
# print samples without Wolbachia
qpcr_no_wolbachia %>%
kbl(booktable=TRUE) %>%
kable_paper("hover", full_width = F)| Sample | Group | |
|---|---|---|
| CTC1 | CTC1 | 0 |
| CTC10 | CTC10 | 0 |
| CTC11 | CTC11 | 0 |
| CTC13 | CTC13 | 0 |
| CTC14 | CTC14 | 0 |
| CTC15 | CTC15 | 0 |
| CTC2 | CTC2 | 0 |
| CTC3 | CTC3 | 0 |
| CTC4 | CTC4 | 0 |
| CTC6 | CTC6 | 0 |
| CTC7 | CTC7 | 0 |
# print samples with Wolbachia
qpcr_wolbachia %>%
kbl(booktable=TRUE) %>%
kable_paper("hover", full_width = F)| Sample | Ratio_Target_Ref | |
|---|---|---|
| CTC12 | CTC12 | 0.000197 |
| NJ39 | NJ39 | 1.395000 |
| NJ41 | NJ41 | 3.481000 |
| NJ42 | NJ42 | 1.912000 |
| NJ46 | NJ46 | 0.496800 |
| S18 | S18 | 5.689000 |
| S19 | S19 | 13.310000 |
| S20 | S20 | 9.815000 |
| S21 | S21 | 4.573000 |
| S22 | S22 | 2.879000 |
| S23 | S23 | 7.873000 |
| S25 | S25 | 4.360000 |
| S27 | S27 | 0.114000 |
| S30 | S30 | 5.614000 |
| S42 | S42 | 2.207000 |
| S43 | S43 | 0.017600 |
| S45 | S45 | 0.352100 |
| S47 | S47 | 0.312800 |
| S48 | S48 | 0.237800 |
| S50 | S50 | 0.209600 |
| S51 | S51 | 3.615000 |
| S52 | S52 | 1.272000 |
| S64 | S64 | 2.427000 |
# add column with true ID in qpcr table
setwd(path_qpcr)
metadata <- openxlsx::read.xlsx("Supplementary_Table_1_V4_Shorten_w_ReadDepth_CheckqPCR.xlsx", sheet=1)
qpcr$True_TubeID <- rownames(qpcr)
qpcr <- qpcr %>% merge(metadata, by="True_TubeID", all.x=TRUE)
qpcr <- qpcr %>% select(c(True_TubeID, Sample.ID, Ratio_Target_Ref))
# print the qpcr table
qpcr %>%
kbl(booktable=TRUE) %>%
kable_paper("hover", full_width = F)| True_TubeID | Sample.ID | Ratio_Target_Ref |
|---|---|---|
| CTC1 | CTC1 | 0.000000 |
| CTC10 | CTC10 | 0.000000 |
| CTC11 | CTC11 | 0.000000 |
| CTC12 | CTC12 | 0.000197 |
| CTC13 | CTC13 | 0.000000 |
| CTC14 | CTC14 | 0.000000 |
| CTC15 | CTC15 | 0.000000 |
| CTC2 | CTC2 | 0.000000 |
| CTC3 | CTC3 | 0.000000 |
| CTC4 | CTC4 | 0.000000 |
| CTC6 | CTC6 | 0.000000 |
| CTC7 | CTC7 | 0.000000 |
| NJ39 | NP27 | 1.395000 |
| NJ41 | NP29 | 3.481000 |
| NJ42 | NP30 | 1.912000 |
| NJ46 | NP34 | 0.496800 |
| S18 | S18 | 5.689000 |
| S19 | S19 | 13.310000 |
| S20 | S20 | 9.815000 |
| S21 | S21 | 4.573000 |
| S22 | S22 | 2.879000 |
| S23 | S23 | 7.873000 |
| S25 | S25 | 4.360000 |
| S27 | S27 | 0.114000 |
| S30 | S30 | 5.614000 |
| S42 | S42 | 2.207000 |
| S43 | S43 | 0.017600 |
| S45 | S45 | 0.352100 |
| S47 | S47 | 0.312800 |
| S48 | S48 | 0.237800 |
| S50 | S50 | 0.209600 |
| S51 | S51 | 3.615000 |
| S52 | S52 | 1.272000 |
| S64 | S64 | 2.427000 |
HCA
Run HCA
Dendogram
Make a dataframe from these groups
Create dataframe
## CTC12 NJ39 NJ41 NJ42 NJ46 S18 S21 S22 S25 S27 S30 S42 S43
## 1 1 1 1 1 1 1 1 1 1 1 1 1
## S45 S47 S48 S50 S51 S52 S64 S19 S20 S23
## 1 1 1 1 1 1 1 2 2 2
Add samples with no Wolbachia
qpcr.groups <- qpcr.groups %>% rbind(qpcr_no_wolbachia)
# print
qpcr.groups %>%
kbl(booktable=TRUE) %>%
kable_paper("hover", full_width = F)| Group | Sample | |
|---|---|---|
| CTC12 | 1 | CTC12 |
| NJ39 | 1 | NJ39 |
| NJ41 | 1 | NJ41 |
| NJ42 | 1 | NJ42 |
| NJ46 | 1 | NJ46 |
| S18 | 1 | S18 |
| S21 | 1 | S21 |
| S22 | 1 | S22 |
| S25 | 1 | S25 |
| S27 | 1 | S27 |
| S30 | 1 | S30 |
| S42 | 1 | S42 |
| S43 | 1 | S43 |
| S45 | 1 | S45 |
| S47 | 1 | S47 |
| S48 | 1 | S48 |
| S50 | 1 | S50 |
| S51 | 1 | S51 |
| S52 | 1 | S52 |
| S64 | 1 | S64 |
| S19 | 2 | S19 |
| S20 | 2 | S20 |
| S23 | 2 | S23 |
| CTC1 | 0 | CTC1 |
| CTC10 | 0 | CTC10 |
| CTC11 | 0 | CTC11 |
| CTC13 | 0 | CTC13 |
| CTC14 | 0 | CTC14 |
| CTC15 | 0 | CTC15 |
| CTC2 | 0 | CTC2 |
| CTC3 | 0 | CTC3 |
| CTC4 | 0 | CTC4 |
| CTC6 | 0 | CTC6 |
| CTC7 | 0 | CTC7 |
Add groups in the qpcr table
rownames(qpcr) <- qpcr$True_TubeID
qpcr <- qpcr %>% merge(qpcr.groups, by="row.names")
rownames(qpcr) <- qpcr$Sample.ID
qpcr <- qpcr[,-1]
qpcr <- qpcr %>% select(-c(Sample))
# print the qpcr table
qpcr %>%
kbl(booktable=TRUE) %>%
kable_paper("hover", full_width = F)| True_TubeID | Sample.ID | Ratio_Target_Ref | Group | |
|---|---|---|---|---|
| CTC1 | CTC1 | CTC1 | 0.000000 | 0 |
| CTC10 | CTC10 | CTC10 | 0.000000 | 0 |
| CTC11 | CTC11 | CTC11 | 0.000000 | 0 |
| CTC12 | CTC12 | CTC12 | 0.000197 | 1 |
| CTC13 | CTC13 | CTC13 | 0.000000 | 0 |
| CTC14 | CTC14 | CTC14 | 0.000000 | 0 |
| CTC15 | CTC15 | CTC15 | 0.000000 | 0 |
| CTC2 | CTC2 | CTC2 | 0.000000 | 0 |
| CTC3 | CTC3 | CTC3 | 0.000000 | 0 |
| CTC4 | CTC4 | CTC4 | 0.000000 | 0 |
| CTC6 | CTC6 | CTC6 | 0.000000 | 0 |
| CTC7 | CTC7 | CTC7 | 0.000000 | 0 |
| NP27 | NJ39 | NP27 | 1.395000 | 1 |
| NP29 | NJ41 | NP29 | 3.481000 | 1 |
| NP30 | NJ42 | NP30 | 1.912000 | 1 |
| NP34 | NJ46 | NP34 | 0.496800 | 1 |
| S18 | S18 | S18 | 5.689000 | 1 |
| S19 | S19 | S19 | 13.310000 | 2 |
| S20 | S20 | S20 | 9.815000 | 2 |
| S21 | S21 | S21 | 4.573000 | 1 |
| S22 | S22 | S22 | 2.879000 | 1 |
| S23 | S23 | S23 | 7.873000 | 2 |
| S25 | S25 | S25 | 4.360000 | 1 |
| S27 | S27 | S27 | 0.114000 | 1 |
| S30 | S30 | S30 | 5.614000 | 1 |
| S42 | S42 | S42 | 2.207000 | 1 |
| S43 | S43 | S43 | 0.017600 | 1 |
| S45 | S45 | S45 | 0.352100 | 1 |
| S47 | S47 | S47 | 0.312800 | 1 |
| S48 | S48 | S48 | 0.237800 | 1 |
| S50 | S50 | S50 | 0.209600 | 1 |
| S51 | S51 | S51 | 3.615000 | 1 |
| S52 | S52 | S52 | 1.272000 | 1 |
| S64 | S64 | S64 | 2.427000 | 1 |
Groups statistics
With raw labels
# make table with stats
qpcr.summary <- summary.table(qpcr, 0, "Ratio_Target_Ref") %>% cbind(summary.table(qpcr, 1, "Ratio_Target_Ref")) %>% cbind(summary.table(qpcr, 2, "Ratio_Target_Ref"))
# print
qpcr.summary %>%
kbl(booktable=TRUE) %>%
kable_paper("hover", full_width = F)| 0 | 1 | 2 | |
|---|---|---|---|
| Min. | 0 | 0.000197 | 7.87300 |
| 1st Qu. | 0 | 0.294050 | 8.84400 |
| Median | 0 | 1.653500 | 9.81500 |
| Mean | 0 | 2.058245 | 10.33267 |
| 3rd Qu. | 0 | 3.514500 | 11.56250 |
| Max. | 0 | 5.689000 | 13.31000 |
With new labels
## [1] "0" "1" "2"
# change levels of Group
qpcr[qpcr$Group==0, "Group"] <- "No infection"
qpcr[qpcr$Group==1, "Group"] <- "Low infection"
qpcr[qpcr$Group==2, "Group"] <- "High infection"
# make table with stats
qpcr.summary.labels <- summary.table(qpcr, "No infection", "Ratio_Target_Ref") %>% cbind(summary.table(qpcr, "Low infection", "Ratio_Target_Ref")) %>% cbind(summary.table(qpcr, "High infection", "Ratio_Target_Ref"))
# print
qpcr.summary.labels %>%
kbl(booktable=TRUE) %>%
kable_paper("hover", full_width = F)| No infection | Low infection | High infection | |
|---|---|---|---|
| Min. | 0 | 0.000197 | 7.87300 |
| 1st Qu. | 0 | 0.294050 | 8.84400 |
| Median | 0 | 1.653500 | 9.81500 |
| Mean | 0 | 2.058245 | 10.33267 |
| 3rd Qu. | 0 | 3.514500 | 11.56250 |
| Max. | 0 | 5.689000 | 13.31000 |
NMDS (Non-metric MultiDimensional Scaling)
Preparation
Merge groups of data from phyloseq and data from qpcr
# print wolbachia.groups (phyloseq)
wolbachia.groups %>%
kbl(booktable=TRUE) %>%
kable_paper("hover", full_width = F)| Sample | Group | Wolbachia | |
|---|---|---|---|
| CTC1 | CTC1 | Low infection | 0.0000000 |
| CTC10 | CTC10 | Low infection | 0.0352626 |
| CTC11 | CTC11 | Low infection | 0.0077123 |
| CTC12 | CTC12 | Low infection | 0.0191972 |
| CTC13 | CTC13 | Low infection | 0.0026615 |
| CTC14 | CTC14 | Medium infection | 0.1074139 |
| CTC15 | CTC15 | Low infection | 0.0094047 |
| CTC2 | CTC2 | Low infection | 0.0087316 |
| CTC3 | CTC3 | Low infection | 0.0007014 |
| CTC4 | CTC4 | Low infection | 0.0088827 |
| CTC5 | CTC5 | Low infection | 0.0017734 |
| CTC6 | CTC6 | Low infection | 0.0020249 |
| CTC7 | CTC7 | Low infection | 0.0050378 |
| CTC9 | CTC9 | Low infection | 0.0058233 |
| NP14 | NP14 | High infection | 0.9888611 |
| NP2 | NP2 | High infection | 0.9993923 |
| NP20 | NP20 | Medium infection | 0.4705882 |
| NP27 | NP27 | High infection | 0.9789303 |
| NP29 | NP29 | Low infection | 0.0492611 |
| NP30 | NP30 | Medium infection | 0.2894737 |
| NP34 | NP34 | Medium infection | 0.3157895 |
| NP35 | NP35 | Low infection | 0.0032998 |
| NP36 | NP36 | Low infection | 0.0000000 |
| NP37 | NP37 | Low infection | 0.0002480 |
| NP38 | NP38 | Low infection | 0.0218636 |
| NP39 | NP39 | Low infection | 0.0023176 |
| NP41 | NP41 | Low infection | 0.0000162 |
| NP42 | NP42 | Low infection | 0.0004591 |
| NP43 | NP43 | Low infection | 0.0070779 |
| NP44 | NP44 | Low infection | 0.0000127 |
| NP5 | NP5 | High infection | 0.9999973 |
| NP8 | NP8 | High infection | 0.9999492 |
| S100 | S100 | High infection | 1.0000000 |
| S101 | S101 | High infection | 0.8433759 |
| S102 | S102 | High infection | 1.0000000 |
| S104 | S104 | High infection | 0.9999237 |
| S105 | S105 | High infection | 0.9999100 |
| S106 | S106 | High infection | 1.0000000 |
| S107 | S107 | High infection | 0.9999808 |
| S108 | S108 | High infection | 0.9999641 |
| S109 | S109 | High infection | 0.9999492 |
| S110 | S110 | High infection | 1.0000000 |
| S113 | S113 | High infection | 0.9941796 |
| S121 | S121 | High infection | 1.0000000 |
| S122 | S122 | High infection | 1.0000000 |
| S123 | S123 | High infection | 1.0000000 |
| S124 | S124 | High infection | 0.9999449 |
| S126 | S126 | High infection | 0.9997374 |
| S127 | S127 | High infection | 1.0000000 |
| S128 | S128 | High infection | 1.0000000 |
| S146 | S146 | High infection | 0.9993204 |
| S147 | S147 | High infection | 0.9999205 |
| S148 | S148 | High infection | 0.9993646 |
| S150 | S150 | High infection | 1.0000000 |
| S151 | S151 | High infection | 1.0000000 |
| S152 | S152 | High infection | 0.9998674 |
| S153 | S153 | High infection | 0.9998826 |
| S154 | S154 | High infection | 0.9993767 |
| S160 | S160 | High infection | 0.9999035 |
| S162 | S162 | High infection | 0.9999603 |
| S163 | S163 | High infection | 0.9999188 |
| S164 | S164 | High infection | 0.9993741 |
| S165 | S165 | High infection | 0.9995488 |
| S166 | S166 | High infection | 0.9992847 |
| S167 | S167 | High infection | 0.9997153 |
| S169 | S169 | High infection | 1.0000000 |
| S170 | S170 | High infection | 1.0000000 |
| S18 | S18 | High infection | 0.9909091 |
| S19 | S19 | Medium infection | 0.6928881 |
| S20 | S20 | High infection | 0.9376633 |
| S21 | S21 | Medium infection | 0.6642608 |
| S22 | S22 | Medium infection | 0.5062214 |
| S23 | S23 | High infection | 0.9421566 |
| S24 | S24 | High infection | 0.9054080 |
| S25 | S25 | High infection | 0.9970647 |
| S26 | S26 | High infection | 0.9733037 |
| S27 | S27 | Medium infection | 0.2155884 |
| S28 | S28 | High infection | 0.9646902 |
| S29 | S29 | Low infection | 0.0098804 |
| S30 | S30 | Medium infection | 0.1662131 |
| S31 | S31 | Medium infection | 0.1847797 |
| S32 | S32 | Medium infection | 0.2115416 |
| S33 | S33 | Medium infection | 0.3796025 |
| S34 | S34 | Medium infection | 0.5882353 |
| S35 | S35 | Medium infection | 0.4507416 |
| S36 | S36 | Medium infection | 0.1426408 |
| S37 | S37 | Low infection | 0.0420938 |
| S38 | S38 | High infection | 0.8371053 |
| S39 | S39 | Medium infection | 0.3274431 |
| S40 | S40 | Medium infection | 0.4336303 |
| S41 | S41 | Medium infection | 0.3043478 |
| S42 | S42 | High infection | 0.9997565 |
| S43 | S43 | High infection | 0.9369888 |
| S44 | S44 | High infection | 0.9795664 |
| S45 | S45 | High infection | 0.9815447 |
| S46 | S46 | Medium infection | 0.4921260 |
| S47 | S47 | Medium infection | 0.4081812 |
| S48 | S48 | High infection | 0.9938841 |
| S49 | S49 | High infection | 0.9925667 |
| S50 | S50 | High infection | 0.9898880 |
| S51 | S51 | High infection | 0.9936881 |
| S52 | S52 | High infection | 0.9833898 |
| S53 | S53 | High infection | 1.0000000 |
| S55 | S55 | High infection | 0.9999604 |
| S56 | S56 | High infection | 1.0000000 |
| S57 | S57 | High infection | 0.9998373 |
| S58 | S58 | High infection | 0.9977978 |
| S59 | S59 | High infection | 1.0000000 |
| S60 | S60 | High infection | 1.0000000 |
| S61 | S61 | High infection | 1.0000000 |
| S63 | S63 | High infection | 0.9955108 |
| S64 | S64 | High infection | 1.0000000 |
| S65 | S65 | Medium infection | 0.3235294 |
| S77 | S77 | Low infection | 0.0000000 |
| S79 | S79 | High infection | 1.0000000 |
| S80 | S80 | High infection | 1.0000000 |
| S83 | S83 | High infection | 1.0000000 |
| S84 | S84 | High infection | 1.0000000 |
| S85 | S85 | High infection | 0.9994486 |
| S86 | S86 | High infection | 1.0000000 |
| S87 | S87 | High infection | 0.9998787 |
| S88 | S88 | High infection | 0.9977787 |
| S89 | S89 | Low infection | 0.0032180 |
| True_TubeID | Sample.ID | Ratio_Target_Ref | Group | |
|---|---|---|---|---|
| CTC1 | CTC1 | CTC1 | 0.000000 | No infection |
| CTC10 | CTC10 | CTC10 | 0.000000 | No infection |
| CTC11 | CTC11 | CTC11 | 0.000000 | No infection |
| CTC12 | CTC12 | CTC12 | 0.000197 | Low infection |
| CTC13 | CTC13 | CTC13 | 0.000000 | No infection |
| CTC14 | CTC14 | CTC14 | 0.000000 | No infection |
| CTC15 | CTC15 | CTC15 | 0.000000 | No infection |
| CTC2 | CTC2 | CTC2 | 0.000000 | No infection |
| CTC3 | CTC3 | CTC3 | 0.000000 | No infection |
| CTC4 | CTC4 | CTC4 | 0.000000 | No infection |
| CTC6 | CTC6 | CTC6 | 0.000000 | No infection |
| CTC7 | CTC7 | CTC7 | 0.000000 | No infection |
| NP27 | NJ39 | NP27 | 1.395000 | Low infection |
| NP29 | NJ41 | NP29 | 3.481000 | Low infection |
| NP30 | NJ42 | NP30 | 1.912000 | Low infection |
| NP34 | NJ46 | NP34 | 0.496800 | Low infection |
| S18 | S18 | S18 | 5.689000 | Low infection |
| S19 | S19 | S19 | 13.310000 | High infection |
| S20 | S20 | S20 | 9.815000 | High infection |
| S21 | S21 | S21 | 4.573000 | Low infection |
| S22 | S22 | S22 | 2.879000 | Low infection |
| S23 | S23 | S23 | 7.873000 | High infection |
| S25 | S25 | S25 | 4.360000 | Low infection |
| S27 | S27 | S27 | 0.114000 | Low infection |
| S30 | S30 | S30 | 5.614000 | Low infection |
| S42 | S42 | S42 | 2.207000 | Low infection |
| S43 | S43 | S43 | 0.017600 | Low infection |
| S45 | S45 | S45 | 0.352100 | Low infection |
| S47 | S47 | S47 | 0.312800 | Low infection |
| S48 | S48 | S48 | 0.237800 | Low infection |
| S50 | S50 | S50 | 0.209600 | Low infection |
| S51 | S51 | S51 | 3.615000 | Low infection |
| S52 | S52 | S52 | 1.272000 | Low infection |
| S64 | S64 | S64 | 2.427000 | Low infection |
# merge
colnames(wolbachia.groups) <- c("Sample", "Infection", "Wolbachia")
colnames(qpcr) <- c("qPCR_ID", "Sample", "qPCR_Ratio_Target_Ref", "qPCR_Infection")
all.groups <- wolbachia.groups %>% merge(qpcr, by="Sample", all.x=TRUE)
all.groups %>%
kbl(booktable=TRUE) %>%
kable_paper("hover", full_width = F)| Sample | Infection | Wolbachia | qPCR_ID | qPCR_Ratio_Target_Ref | qPCR_Infection |
|---|---|---|---|---|---|
| CTC1 | Low infection | 0.0000000 | CTC1 | 0.000000 | No infection |
| CTC10 | Low infection | 0.0352626 | CTC10 | 0.000000 | No infection |
| CTC11 | Low infection | 0.0077123 | CTC11 | 0.000000 | No infection |
| CTC12 | Low infection | 0.0191972 | CTC12 | 0.000197 | Low infection |
| CTC13 | Low infection | 0.0026615 | CTC13 | 0.000000 | No infection |
| CTC14 | Medium infection | 0.1074139 | CTC14 | 0.000000 | No infection |
| CTC15 | Low infection | 0.0094047 | CTC15 | 0.000000 | No infection |
| CTC2 | Low infection | 0.0087316 | CTC2 | 0.000000 | No infection |
| CTC3 | Low infection | 0.0007014 | CTC3 | 0.000000 | No infection |
| CTC4 | Low infection | 0.0088827 | CTC4 | 0.000000 | No infection |
| CTC5 | Low infection | 0.0017734 | NA | NA | NA |
| CTC6 | Low infection | 0.0020249 | CTC6 | 0.000000 | No infection |
| CTC7 | Low infection | 0.0050378 | CTC7 | 0.000000 | No infection |
| CTC9 | Low infection | 0.0058233 | NA | NA | NA |
| NP14 | High infection | 0.9888611 | NA | NA | NA |
| NP2 | High infection | 0.9993923 | NA | NA | NA |
| NP20 | Medium infection | 0.4705882 | NA | NA | NA |
| NP27 | High infection | 0.9789303 | NJ39 | 1.395000 | Low infection |
| NP29 | Low infection | 0.0492611 | NJ41 | 3.481000 | Low infection |
| NP30 | Medium infection | 0.2894737 | NJ42 | 1.912000 | Low infection |
| NP34 | Medium infection | 0.3157895 | NJ46 | 0.496800 | Low infection |
| NP35 | Low infection | 0.0032998 | NA | NA | NA |
| NP36 | Low infection | 0.0000000 | NA | NA | NA |
| NP37 | Low infection | 0.0002480 | NA | NA | NA |
| NP38 | Low infection | 0.0218636 | NA | NA | NA |
| NP39 | Low infection | 0.0023176 | NA | NA | NA |
| NP41 | Low infection | 0.0000162 | NA | NA | NA |
| NP42 | Low infection | 0.0004591 | NA | NA | NA |
| NP43 | Low infection | 0.0070779 | NA | NA | NA |
| NP44 | Low infection | 0.0000127 | NA | NA | NA |
| NP5 | High infection | 0.9999973 | NA | NA | NA |
| NP8 | High infection | 0.9999492 | NA | NA | NA |
| S100 | High infection | 1.0000000 | NA | NA | NA |
| S101 | High infection | 0.8433759 | NA | NA | NA |
| S102 | High infection | 1.0000000 | NA | NA | NA |
| S104 | High infection | 0.9999237 | NA | NA | NA |
| S105 | High infection | 0.9999100 | NA | NA | NA |
| S106 | High infection | 1.0000000 | NA | NA | NA |
| S107 | High infection | 0.9999808 | NA | NA | NA |
| S108 | High infection | 0.9999641 | NA | NA | NA |
| S109 | High infection | 0.9999492 | NA | NA | NA |
| S110 | High infection | 1.0000000 | NA | NA | NA |
| S113 | High infection | 0.9941796 | NA | NA | NA |
| S121 | High infection | 1.0000000 | NA | NA | NA |
| S122 | High infection | 1.0000000 | NA | NA | NA |
| S123 | High infection | 1.0000000 | NA | NA | NA |
| S124 | High infection | 0.9999449 | NA | NA | NA |
| S126 | High infection | 0.9997374 | NA | NA | NA |
| S127 | High infection | 1.0000000 | NA | NA | NA |
| S128 | High infection | 1.0000000 | NA | NA | NA |
| S146 | High infection | 0.9993204 | NA | NA | NA |
| S147 | High infection | 0.9999205 | NA | NA | NA |
| S148 | High infection | 0.9993646 | NA | NA | NA |
| S150 | High infection | 1.0000000 | NA | NA | NA |
| S151 | High infection | 1.0000000 | NA | NA | NA |
| S152 | High infection | 0.9998674 | NA | NA | NA |
| S153 | High infection | 0.9998826 | NA | NA | NA |
| S154 | High infection | 0.9993767 | NA | NA | NA |
| S160 | High infection | 0.9999035 | NA | NA | NA |
| S162 | High infection | 0.9999603 | NA | NA | NA |
| S163 | High infection | 0.9999188 | NA | NA | NA |
| S164 | High infection | 0.9993741 | NA | NA | NA |
| S165 | High infection | 0.9995488 | NA | NA | NA |
| S166 | High infection | 0.9992847 | NA | NA | NA |
| S167 | High infection | 0.9997153 | NA | NA | NA |
| S169 | High infection | 1.0000000 | NA | NA | NA |
| S170 | High infection | 1.0000000 | NA | NA | NA |
| S18 | High infection | 0.9909091 | S18 | 5.689000 | Low infection |
| S19 | Medium infection | 0.6928881 | S19 | 13.310000 | High infection |
| S20 | High infection | 0.9376633 | S20 | 9.815000 | High infection |
| S21 | Medium infection | 0.6642608 | S21 | 4.573000 | Low infection |
| S22 | Medium infection | 0.5062214 | S22 | 2.879000 | Low infection |
| S23 | High infection | 0.9421566 | S23 | 7.873000 | High infection |
| S24 | High infection | 0.9054080 | NA | NA | NA |
| S25 | High infection | 0.9970647 | S25 | 4.360000 | Low infection |
| S26 | High infection | 0.9733037 | NA | NA | NA |
| S27 | Medium infection | 0.2155884 | S27 | 0.114000 | Low infection |
| S28 | High infection | 0.9646902 | NA | NA | NA |
| S29 | Low infection | 0.0098804 | NA | NA | NA |
| S30 | Medium infection | 0.1662131 | S30 | 5.614000 | Low infection |
| S31 | Medium infection | 0.1847797 | NA | NA | NA |
| S32 | Medium infection | 0.2115416 | NA | NA | NA |
| S33 | Medium infection | 0.3796025 | NA | NA | NA |
| S34 | Medium infection | 0.5882353 | NA | NA | NA |
| S35 | Medium infection | 0.4507416 | NA | NA | NA |
| S36 | Medium infection | 0.1426408 | NA | NA | NA |
| S37 | Low infection | 0.0420938 | NA | NA | NA |
| S38 | High infection | 0.8371053 | NA | NA | NA |
| S39 | Medium infection | 0.3274431 | NA | NA | NA |
| S40 | Medium infection | 0.4336303 | NA | NA | NA |
| S41 | Medium infection | 0.3043478 | NA | NA | NA |
| S42 | High infection | 0.9997565 | S42 | 2.207000 | Low infection |
| S43 | High infection | 0.9369888 | S43 | 0.017600 | Low infection |
| S44 | High infection | 0.9795664 | NA | NA | NA |
| S45 | High infection | 0.9815447 | S45 | 0.352100 | Low infection |
| S46 | Medium infection | 0.4921260 | NA | NA | NA |
| S47 | Medium infection | 0.4081812 | S47 | 0.312800 | Low infection |
| S48 | High infection | 0.9938841 | S48 | 0.237800 | Low infection |
| S49 | High infection | 0.9925667 | NA | NA | NA |
| S50 | High infection | 0.9898880 | S50 | 0.209600 | Low infection |
| S51 | High infection | 0.9936881 | S51 | 3.615000 | Low infection |
| S52 | High infection | 0.9833898 | S52 | 1.272000 | Low infection |
| S53 | High infection | 1.0000000 | NA | NA | NA |
| S55 | High infection | 0.9999604 | NA | NA | NA |
| S56 | High infection | 1.0000000 | NA | NA | NA |
| S57 | High infection | 0.9998373 | NA | NA | NA |
| S58 | High infection | 0.9977978 | NA | NA | NA |
| S59 | High infection | 1.0000000 | NA | NA | NA |
| S60 | High infection | 1.0000000 | NA | NA | NA |
| S61 | High infection | 1.0000000 | NA | NA | NA |
| S63 | High infection | 0.9955108 | NA | NA | NA |
| S64 | High infection | 1.0000000 | S64 | 2.427000 | Low infection |
| S65 | Medium infection | 0.3235294 | NA | NA | NA |
| S77 | Low infection | 0.0000000 | NA | NA | NA |
| S79 | High infection | 1.0000000 | NA | NA | NA |
| S80 | High infection | 1.0000000 | NA | NA | NA |
| S83 | High infection | 1.0000000 | NA | NA | NA |
| S84 | High infection | 1.0000000 | NA | NA | NA |
| S85 | High infection | 0.9994486 | NA | NA | NA |
| S86 | High infection | 1.0000000 | NA | NA | NA |
| S87 | High infection | 0.9998787 | NA | NA | NA |
| S88 | High infection | 0.9977787 | NA | NA | NA |
| S89 | Low infection | 0.0032180 | NA | NA | NA |
# optionnal : all.groups for comparison
all.groups.comparison <- wolbachia.groups %>% merge(qpcr, by="Sample")
all.groups.comparison <- all.groups.comparison %>% select(c(Sample, qPCR_ID, Infection, Wolbachia, qPCR_Ratio_Target_Ref, qPCR_Infection))
all.groups.comparison$Infection <- all.groups.comparison$Infection %>% as.factor()
all.groups.comparison$Infection <- factor(all.groups.comparison$Infection, levels = c("Low infection", "Medium infection", "High infection"))
all.groups.comparison <- all.groups.comparison[order(all.groups.comparison$Infection),]
all.groups.comparison %>%
kbl(booktable=TRUE) %>%
kable_paper("hover", full_width = F)| Sample | qPCR_ID | Infection | Wolbachia | qPCR_Ratio_Target_Ref | qPCR_Infection | |
|---|---|---|---|---|---|---|
| 1 | CTC1 | CTC1 | Low infection | 0.0000000 | 0.000000 | No infection |
| 2 | CTC10 | CTC10 | Low infection | 0.0352626 | 0.000000 | No infection |
| 3 | CTC11 | CTC11 | Low infection | 0.0077123 | 0.000000 | No infection |
| 4 | CTC12 | CTC12 | Low infection | 0.0191972 | 0.000197 | Low infection |
| 5 | CTC13 | CTC13 | Low infection | 0.0026615 | 0.000000 | No infection |
| 7 | CTC15 | CTC15 | Low infection | 0.0094047 | 0.000000 | No infection |
| 8 | CTC2 | CTC2 | Low infection | 0.0087316 | 0.000000 | No infection |
| 9 | CTC3 | CTC3 | Low infection | 0.0007014 | 0.000000 | No infection |
| 10 | CTC4 | CTC4 | Low infection | 0.0088827 | 0.000000 | No infection |
| 11 | CTC6 | CTC6 | Low infection | 0.0020249 | 0.000000 | No infection |
| 12 | CTC7 | CTC7 | Low infection | 0.0050378 | 0.000000 | No infection |
| 14 | NP29 | NJ41 | Low infection | 0.0492611 | 3.481000 | Low infection |
| 6 | CTC14 | CTC14 | Medium infection | 0.1074139 | 0.000000 | No infection |
| 15 | NP30 | NJ42 | Medium infection | 0.2894737 | 1.912000 | Low infection |
| 16 | NP34 | NJ46 | Medium infection | 0.3157895 | 0.496800 | Low infection |
| 18 | S19 | S19 | Medium infection | 0.6928881 | 13.310000 | High infection |
| 20 | S21 | S21 | Medium infection | 0.6642608 | 4.573000 | Low infection |
| 21 | S22 | S22 | Medium infection | 0.5062214 | 2.879000 | Low infection |
| 24 | S27 | S27 | Medium infection | 0.2155884 | 0.114000 | Low infection |
| 25 | S30 | S30 | Medium infection | 0.1662131 | 5.614000 | Low infection |
| 29 | S47 | S47 | Medium infection | 0.4081812 | 0.312800 | Low infection |
| 13 | NP27 | NJ39 | High infection | 0.9789303 | 1.395000 | Low infection |
| 17 | S18 | S18 | High infection | 0.9909091 | 5.689000 | Low infection |
| 19 | S20 | S20 | High infection | 0.9376633 | 9.815000 | High infection |
| 22 | S23 | S23 | High infection | 0.9421566 | 7.873000 | High infection |
| 23 | S25 | S25 | High infection | 0.9970647 | 4.360000 | Low infection |
| 26 | S42 | S42 | High infection | 0.9997565 | 2.207000 | Low infection |
| 27 | S43 | S43 | High infection | 0.9369888 | 0.017600 | Low infection |
| 28 | S45 | S45 | High infection | 0.9815447 | 0.352100 | Low infection |
| 30 | S48 | S48 | High infection | 0.9938841 | 0.237800 | Low infection |
| 31 | S50 | S50 | High infection | 0.9898880 | 0.209600 | Low infection |
| 32 | S51 | S51 | High infection | 0.9936881 | 3.615000 | Low infection |
| 33 | S52 | S52 | High infection | 0.9833898 | 1.272000 | Low infection |
| 34 | S64 | S64 | High infection | 1.0000000 | 2.427000 | Low infection |
all.groups.comparison <- all.groups.comparison %>% merge(sam %>% select(c(Sample, Strain, Species)), by="Sample")
all.groups.comparison <- all.groups.comparison[order(all.groups.comparison$Infection),]
all.groups.comparison %>%
kbl(booktable=TRUE) %>%
kable_paper("hover", full_width = F)| Sample | qPCR_ID | Infection | Wolbachia | qPCR_Ratio_Target_Ref | qPCR_Infection | Strain | Species | |
|---|---|---|---|---|---|---|---|---|
| 1 | CTC1 | CTC1 | Low infection | 0.0000000 | 0.000000 | No infection | Laboratory - Slab TC (Wolbachia -) | Culex quinquefasciatus |
| 2 | CTC10 | CTC10 | Low infection | 0.0352626 | 0.000000 | No infection | Laboratory - Slab TC (Wolbachia -) | Culex quinquefasciatus |
| 3 | CTC11 | CTC11 | Low infection | 0.0077123 | 0.000000 | No infection | Laboratory - Slab TC (Wolbachia -) | Culex quinquefasciatus |
| 4 | CTC12 | CTC12 | Low infection | 0.0191972 | 0.000197 | Low infection | Laboratory - Slab TC (Wolbachia -) | Culex quinquefasciatus |
| 5 | CTC13 | CTC13 | Low infection | 0.0026615 | 0.000000 | No infection | Laboratory - Slab TC (Wolbachia -) | Culex quinquefasciatus |
| 7 | CTC15 | CTC15 | Low infection | 0.0094047 | 0.000000 | No infection | Laboratory - Slab TC (Wolbachia -) | Culex quinquefasciatus |
| 8 | CTC2 | CTC2 | Low infection | 0.0087316 | 0.000000 | No infection | Laboratory - Slab TC (Wolbachia -) | Culex quinquefasciatus |
| 9 | CTC3 | CTC3 | Low infection | 0.0007014 | 0.000000 | No infection | Laboratory - Slab TC (Wolbachia -) | Culex quinquefasciatus |
| 10 | CTC4 | CTC4 | Low infection | 0.0088827 | 0.000000 | No infection | Laboratory - Slab TC (Wolbachia -) | Culex quinquefasciatus |
| 11 | CTC6 | CTC6 | Low infection | 0.0020249 | 0.000000 | No infection | Laboratory - Slab TC (Wolbachia -) | Culex quinquefasciatus |
| 12 | CTC7 | CTC7 | Low infection | 0.0050378 | 0.000000 | No infection | Laboratory - Slab TC (Wolbachia -) | Culex quinquefasciatus |
| 14 | NP29 | NJ41 | Low infection | 0.0492611 | 3.481000 | Low infection | Field - Guadeloupe | Culex quinquefasciatus |
| 6 | CTC14 | CTC14 | Medium infection | 0.1074139 | 0.000000 | No infection | Laboratory - Slab TC (Wolbachia -) | Culex quinquefasciatus |
| 15 | NP30 | NJ42 | Medium infection | 0.2894737 | 1.912000 | Low infection | Field - Guadeloupe | Culex quinquefasciatus |
| 16 | NP34 | NJ46 | Medium infection | 0.3157895 | 0.496800 | Low infection | Field - Guadeloupe | Culex quinquefasciatus |
| 18 | S19 | S19 | Medium infection | 0.6928881 | 13.310000 | High infection | Laboratory - Lavar | Culex pipiens |
| 20 | S21 | S21 | Medium infection | 0.6642608 | 4.573000 | Low infection | Laboratory - Lavar | Culex pipiens |
| 21 | S22 | S22 | Medium infection | 0.5062214 | 2.879000 | Low infection | Laboratory - Lavar | Culex pipiens |
| 24 | S27 | S27 | Medium infection | 0.2155884 | 0.114000 | Low infection | Laboratory - Lavar | Culex pipiens |
| 25 | S30 | S30 | Medium infection | 0.1662131 | 5.614000 | Low infection | Laboratory - Lavar | Culex pipiens |
| 29 | S47 | S47 | Medium infection | 0.4081812 | 0.312800 | Low infection | Field - Camping Europe | Culex pipiens |
| 13 | NP27 | NJ39 | High infection | 0.9789303 | 1.395000 | Low infection | Field - Guadeloupe | Culex quinquefasciatus |
| 17 | S18 | S18 | High infection | 0.9909091 | 5.689000 | Low infection | Laboratory - Lavar | Culex pipiens |
| 19 | S20 | S20 | High infection | 0.9376633 | 9.815000 | High infection | Laboratory - Lavar | Culex pipiens |
| 22 | S23 | S23 | High infection | 0.9421566 | 7.873000 | High infection | Laboratory - Lavar | Culex pipiens |
| 23 | S25 | S25 | High infection | 0.9970647 | 4.360000 | Low infection | Laboratory - Lavar | Culex pipiens |
| 26 | S42 | S42 | High infection | 0.9997565 | 2.207000 | Low infection | Field - Camping Europe | Culex pipiens |
| 27 | S43 | S43 | High infection | 0.9369888 | 0.017600 | Low infection | Field - Camping Europe | Culex pipiens |
| 28 | S45 | S45 | High infection | 0.9815447 | 0.352100 | Low infection | Field - Camping Europe | Culex pipiens |
| 30 | S48 | S48 | High infection | 0.9938841 | 0.237800 | Low infection | Field - Camping Europe | Culex pipiens |
| 31 | S50 | S50 | High infection | 0.9898880 | 0.209600 | Low infection | Field - Bosc | Culex pipiens |
| 32 | S51 | S51 | High infection | 0.9936881 | 3.615000 | Low infection | Field - Bosc | Culex pipiens |
| 33 | S52 | S52 | High infection | 0.9833898 | 1.272000 | Low infection | Field - Bosc | Culex pipiens |
| 34 | S64 | S64 | High infection | 1.0000000 | 2.427000 | Low infection | Field - Bosc | Culex pipiens |
means <- aggregate(qPCR_Ratio_Target_Ref ~ Infection , all.groups.comparison, median)
means$qPCR_Ratio_Target_Ref <- means$qPCR_Ratio_Target_Ref %>% round(2)
ggqqplot(all.groups.comparison$Wolbachia)##
## Shapiro-Wilk normality test
##
## data: all.groups.comparison$Wolbachia
## W = 0.78016, p-value = 1.084e-05
my_comparisons <- list( c("Low infection", "Medium infection"),
c("Medium infection", "High infection"),
c("Low infection", "High infection") )
p1 <- ggplot(all.groups.comparison, aes(x=Infection, y=qPCR_Ratio_Target_Ref, fill=Infection)) +
geom_boxplot(alpha=0.6, outlier.shape = NA)+
geom_jitter(position=position_jitter(0.5), aes(color=Infection), alpha=0.6) +
theme_bw(base_size = 14)+
scale_fill_manual(values = c("blue", "grey", "red"))+
scale_color_manual(values = c("blue", "grey", "red"))+
labs(x=expression(paste("\n", italic("Wolbachia"), "infection groups (from HCA)")),
y="Ratio Target Ref (qPCR)",
fill=expression(paste("\n", italic("Wolbachia"), "infection groups (from HCA)")))+
guides(colour=FALSE)+
#stat_summary(fun.y = "median", geom = "point", shape = 23, size = 3, fill = "white")+
stat_compare_means(comparisons = my_comparisons,
method="wilcox.test",
label="p.signif",
vjust=1.5,
label.y = c(10, 11, 12))+
stat_compare_means(comparisons = my_comparisons,
method="wilcox.test",
label="p.format",
label.y = c(10, 11, 12))+
geom_text(data = means, aes(label = qPCR_Ratio_Target_Ref, y = qPCR_Ratio_Target_Ref + 0.6))
p1p9 <- ggplot(all.groups.comparison, aes(x=Infection, y=Wolbachia, fill=Infection)) +
geom_boxplot(alpha=0.6, outlier.shape = NA)+
geom_jitter(position=position_jitter(0.5), aes(color=Infection), alpha=0.6) +
theme_bw(base_size = 14)+
scale_fill_manual(values = c("blue", "grey", "red"))+
scale_color_manual(values = c("blue", "grey", "red"))+
labs(x=expression(paste("\n", italic("Wolbachia"), "infection groups (from HCA)")),
y=expression(paste("\n", italic("Wolbachia"), "infection (relative abundance)")),
fill="Infection")+
guides(colour=FALSE)+
scale_y_continuous(breaks = seq(0, 1, by = 0.1))
p9setwd(path_tsv)
write.table(all.groups.comparison, "1Ec_qpcr_comparison.tsv", sep="\t", row.names = FALSE, col.names=TRUE)
setwd(path_plot)
tiff("1Ec_boxplot.tiff", units="in", width=10, height=5, res=300)
p1
dev.off()## quartz_off_screen
## 2
## quartz_off_screen
## 2
Add groups in decontam phyloseq object
# load phyloseq object after decontam
setwd(path_rdata)
ps.filter <- readRDS("1D_MED_phyloseq_decontam.rds")
# extract metadata
sam <- data.frame(ps.filter@sam_data)
# add groups in sam table
sam <- sam %>% merge(all.groups %>% select(c(Sample, Infection, Wolbachia, qPCR_Infection, qPCR_Ratio_Target_Ref)), by="Sample", all.x = TRUE)
rownames(sam)<- sam$Sample
# change levels for plot
# levels(sam$Strain) <- c(levels(sam$Strain), "Laboratory - Lavar", "Laboratory - Slab TC (Wolbachia -)")
# sam$Strain[sam$Strain=="Laboratory - Lavar"] <- "Laboratory - Lavar"
# sam$Strain[sam$Strain=="Wolbachia -"] <- "Laboratory - Slab TC (Wolbachia -)"
# sam$Strain <- droplevels(sam$Strain)
# levels(sam$Strain)
# add a column for qPCR measure
sam$qPCR_measure <- "Yes"
sam[is.na(sam$qPCR_Ratio_Target_Ref), "qPCR_measure"] <- "No"
sam[sam$qPCR_measure=="Yes" & sam$qPCR_Ratio_Target_Ref==0, "qPCR_measure"] <- "= 0 (Wolbachia -)"
sam[sam$qPCR_measure=="Yes" & sam$qPCR_Ratio_Target_Ref>0, "qPCR_measure"] <- "> 0 (Wolbachia +)"
# update sam table of phyloseq object
sample_data(ps.filter) <- sam
# transformation percent
ps_percent <- transform_sample_counts(ps.filter, function(x) x/sum(x)*100)
# select whole
ps_percent_whole <- subset_samples(ps_percent, Organ == "Whole")Compute Bray-Curtis distance
set.seed(35)
bray.culex_whole <- ordinate(ps_percent_whole, method="NMDS", distance="bray", trymax=300)## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1088607
## Run 1 stress 0.1070017
## ... New best solution
## ... Procrustes: rmse 0.02445929 max resid 0.1700959
## Run 2 stress 0.1452511
## Run 3 stress 0.1087863
## Run 4 stress 0.1522622
## Run 5 stress 0.1215459
## Run 6 stress 0.1553544
## Run 7 stress 0.119003
## Run 8 stress 0.1259113
## Run 9 stress 0.1064512
## ... New best solution
## ... Procrustes: rmse 0.006708485 max resid 0.03812992
## Run 10 stress 0.1397327
## Run 11 stress 0.1215459
## Run 12 stress 0.1064512
## ... Procrustes: rmse 1.293125e-05 max resid 8.41373e-05
## ... Similar to previous best
## Run 13 stress 0.1521948
## Run 14 stress 0.1389989
## Run 15 stress 0.1558452
## Run 16 stress 0.1090684
## Run 17 stress 0.1174935
## Run 18 stress 0.1090611
## Run 19 stress 0.1070004
## Run 20 stress 0.1175567
## *** Solution reached
Plot
Raw data
# all info
p1 <- plot_ordination(ps_percent_whole, bray.culex_whole, color="Infection", shape=NA) +
labs(x="NMDS1", y = "NMDS2")+
stat_ellipse(geom = "polygon", level=0.95, alpha = 0.25, aes(fill = Species, color=NA))+
scale_fill_manual(values=c("orange", "yellow","green"), labels = c(expression(italic("Aedes aegypti")), expression(italic("Culex pipiens")), expression(italic("Culex quinquefasciatus"))))+
scale_color_manual(values=c("blue", "grey","red"), breaks=c("Low infection", "Medium infection", "High infection"))+
scale_shape_manual(values=c(19,17,15,4,7), labels=c("Field - Bosc", "Field - Camping Europe", "Field - Guadeloupe", "Laboratory - Lavar", expression(paste("Laboratory - Slab TC (", italic("Wolbachia"), "-)"))))+
theme_gray()+
labs(caption=paste("stress:", bray.culex_whole$stress %>% round(2), "\n", "ellipse level: 0.95"))+
#expand_limits(x=c(-2.5, 1.5), y=c(-2, 2.5))+
guides(fill = guide_legend(order = 1),
colour = guide_legend(order = 2),
shape = guide_legend(order = 3))+
labs(fill="Species", col="Wolbachia Group")+
geom_text(mapping = aes(label = Sample), size = 1.5, vjust = 3)+
theme(legend.text.align = 0)## Warning in plot_ordination(ps_percent_whole, bray.culex_whole, color =
## "Infection", : Shape variable was not found in the available data you
## provided.No shape mapped.
p1$layers[[1]] <- NULL
p1[["guides"]][["shape"]][["title"]] <- "Strain"
p1[["guides"]][["colour"]][["title"]] <- expression(paste(italic("Wolbachia"), "Infection"))
p1+
geom_point(aes(shape = Strain), size = 2)# Wolbachia infection
p2 <- plot_ordination(ps_percent_whole, bray.culex_whole, color="Infection", shape="qPCR_measure") +
scale_color_manual(values=c("blue", "#E4A7F5","red"),
breaks=c("Low infection", "Medium infection", "High infection"))+
scale_fill_manual(expression(paste(italic("Wolbachia"), "Infection")),
values=c("blue", "#E4A7F5","red"),
breaks=c("Low infection", "Medium infection", "High infection"))+
scale_shape_manual("qPCR",
values=c(7,2,19),
labels=c(expression(paste("= 0 (", italic("Wolbachia"), "-)")),
expression(paste("> 0 (", italic("Wolbachia"), "+)")),
"No qPCR"))+
labs(x="NMDS1",
y = "NMDS2",
caption=paste("stress:", bray.culex_whole$stress %>% round(2), "\n", "ellipse level: 0.95"))+
geom_text(mapping = aes(label = Sample), size = 1.5, vjust = 3)+
stat_ellipse(geom = "polygon",
alpha = 0.25,
aes(fill=Infection, group=Infection))+
coord_cartesian(xlim = c(-2.8,2.8), ylim=c(-2.5,2.2))+
#expand_limits(x=c(-2.5, 1.5), y=c(-2, 2.5))+
guides(fill = guide_legend(order = 1),
shape = guide_legend(order = 2),
colour=FALSE)+
theme(legend.text.align = 0)
p2# species ("orange", "yellow", "green")
p3 <- plot_ordination(ps_percent_whole, bray.culex_whole, color="Species") +
geom_point(size = 2) +
scale_fill_manual(values=c("#F57462", "#498CF5","#559E36"),
labels = c(expression(italic("Aedes aegypti")),
expression(italic("Culex pipiens")),
expression(italic("Culex quinquefasciatus"))))+
scale_color_manual(values=c("#F57462", "#498CF5","#559E36"))+
labs(x="NMDS1",
y = "NMDS2",
caption=paste("stress:", bray.culex_whole$stress %>% round(2), "\n", "ellipse level: 0.95"))+
geom_text(mapping = aes(label = Sample), size = 1.5, vjust = 3)+
stat_ellipse(geom = "polygon",
alpha = 0.25,
aes(fill = Species))+
coord_cartesian(xlim = c(-2.8,2.8), ylim=c(-2.5,2.2))+
#expand_limits(x=c(-2.5, 1.5), y=c(-2, 2.5))+
guides(fill = guide_legend(order = 1),
colour = FALSE,
shape = guide_legend(order = 3))+
theme(legend.text.align = 0)
p3# Strain
new_location <- c("Field - Guadeloupe", "Laboratory - Slab TC (Wolbachia -)",
"Field - Bosc", "Field - Camping Europe", "Laboratory - Lavar")
ps_percent_whole@sam_data$Strain <- factor(ps_percent_whole@sam_data$Strain, levels = new_location)
p4 <- plot_ordination(ps_percent_whole, bray.culex_whole, color="Strain") +
geom_point(size = 2) +
scale_fill_manual("Strain",
values=c("#00BF7D", "#5B6BF7", "#A3A500", "#F8766D", "#00B0F6"),
labels = c("Field - Guadeloupe",
expression(paste("Laboratory - Slab TC (", italic("Wolbachia"), "-)")),
"Field - Bosc",
"Field - Camping Europe",
"Laboratory - Lavar"))+
scale_color_manual("Strain",
values=c("#00BF7D", "#5B6BF7", "#A3A500", "#F8766D", "#00B0F6"),
labels = c("Field - Guadeloupe",
expression(paste("Laboratory - Slab TC (", italic("Wolbachia"), "-)")),
"Field - Bosc",
"Field - Camping Europe",
"Laboratory - Lavar"))+
labs(x="NMDS1",
y = "NMDS2",
caption=paste("stress:", bray.culex_whole$stress %>% round(2), "\n", "ellipse level: 0.95"))+
stat_ellipse(geom = "polygon", alpha = 0.25, aes(fill=Strain, group=Strain))+
coord_cartesian(xlim = c(-2.8,2.8), ylim=c(-2.5,2.2))+
#expand_limits(x=c(-2.5, 1.5), y=c(-2, 2.5))+
guides(fill = guide_legend(order = 1),
colour = FALSE,
shape = guide_legend(order = 3))+
geom_text(mapping = aes(label = Sample), size = 1.5, vjust = 3)+
theme(legend.text.align = 0)
p4# field
p5 <- plot_ordination(ps_percent_whole, bray.culex_whole, color="Field") +
geom_point(size = 2) +
scale_fill_manual(values=c("brown", "purple" ))+
scale_color_manual(values=c("brown", "purple"))+
labs(x="NMDS1",
y = "NMDS2",
caption=paste("stress:", bray.culex_whole$stress %>% round(2), "\n", "ellipse level: 0.95"))+
stat_ellipse(geom = "polygon", alpha = 0.25, aes(fill=Field, group=Field))+
#coord_cartesian(xlim = c(-2,2.8), ylim=c(-2,1.5))+
coord_cartesian(xlim = c(-2.8,2.8), ylim=c(-2.5,2.2))+
#expand_limits(x=c(-2.5, 1.5), y=c(-2, 2.5))+
guides(fill = guide_legend(order = 1),
colour = FALSE,
shape = guide_legend(order = 3))+
geom_text(mapping = aes(label = Sample), size = 1.5, vjust = 3)
p5# species + strains
p6 <- plot_ordination(ps_percent_whole, bray.culex_whole, color="Strain", shape=NA) +
labs(x="NMDS1", y = "NMDS2")+
scale_color_manual("Strain",
values=c("#00BF7D", "#5B6BF7", "#A3A500", "#F8766D", "#00B0F6"),
labels = c("Field - Guadeloupe",
expression(paste("Laboratory - Slab TC (", italic("Wolbachia"), "-)")),
"Field - Bosc",
"Field - Camping Europe",
"Laboratory - Lavar"))+
scale_shape_manual(values = c(15,16,17), labels=c(expression(italic("Aedes aegypti")), expression(italic("Culex pipiens")), expression(italic("Culex quinquefasciatus"))))+
theme_gray()+
labs(caption=paste("stress:", bray.culex_whole$stress %>% round(2), "\n", "ellipse level: 0.95"))+
guides(fill = guide_legend(order = 1),
colour = guide_legend(order = 2),
shape = guide_legend(order = 3))+
geom_text(mapping = aes(label = Sample), size = 1.5, vjust = 3)+
theme(legend.text.align = 0)+
geom_point(aes(shape = Species_italic), size = 2)## Warning in plot_ordination(ps_percent_whole, bray.culex_whole, color =
## "Strain", : Shape variable was not found in the available data you provided.No
## shape mapped.
sam2 <- sam[sam$Species=="Culex pipiens",]
means <- aggregate(Wolbachia ~ Strain , sam2, median)
means$Wolbachia <- means$Wolbachia %>% round(2)
ggqqplot(sam2$Wolbachia)##
## Shapiro-Wilk normality test
##
## data: sam2$Wolbachia
## W = 0.53983, p-value = 1.266e-14
my_comparisons <- list( c("Field - Bosc", "Field - Camping Europe"),
c("Field - Camping Europe", "Laboratory - Lavar"),
c("Field - Bosc", "Laboratory - Lavar") )
p7 <- ggplot(sam2, aes(x=Strain, y=Wolbachia, fill=Strain)) +
geom_boxplot(alpha=0.6, outlier.shape = NA)+
geom_jitter(position=position_jitter(0.5), aes(color=Strain), alpha=0.6) +
theme_bw(base_size = 14)+
scale_fill_manual(values = c("blue", "grey", "red"))+
scale_color_manual(values = c("blue", "grey", "red"))+
labs(x="\nStrain",
y=expression(paste("\n", italic("Wolbachia"), "infection (relative abundance)")),
fill="Strain")+
guides(colour=FALSE)+
#stat_summary(fun = "mean", geom = "point", shape = 23, size = 3, fill = "white")+
stat_compare_means(comparisons = my_comparisons,
method="wilcox.test",
label="p.signif",
vjust=1.5,
label.y = c(1.1, 1.2, 1.3))+
stat_compare_means(comparisons = my_comparisons,
method="wilcox.test",
label="p.format",
label.y = c(1.1, 1.2, 1.3))+
geom_text(data = means, aes(label = Wolbachia, y = Wolbachia - 0.05))
p7## Warning in wilcox.test.default(c(0.9999999998, 1.0000000001, 1.0000000001, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(c(0.9999999998, 1.0000000001, 1.0000000001, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(c(1.0000000002, 1.0000000002, 0.99992370072, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(c(0.9999999998, 1.0000000001, 1.0000000001, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(c(0.9999999998, 1.0000000001, 1.0000000001, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(c(1.0000000002, 1.0000000002, 0.99992370072, :
## cannot compute exact p-value with ties
means <- aggregate(qPCR_Ratio_Target_Ref ~ Strain , sam2, median)
means$qPCR_Ratio_Target_Ref<- means$qPCR_Ratio_Target_Ref %>% round(2)
wilcox.test(qPCR_Ratio_Target_Ref ~ Strain, sam2[sam2$Strain!="Field - Camping Europe",])##
## Wilcoxon rank sum test
##
## data: qPCR_Ratio_Target_Ref by Strain
## W = 5, p-value = 0.05035
## alternative hypothesis: true location shift is not equal to 0
##
## Shapiro-Wilk normality test
##
## data: sam3$qPCR_Ratio_Target_Ref
## W = 0.86499, p-value = 0.01473
p8 <- ggplot(na.omit(sam2), aes(x=Strain, y=qPCR_Ratio_Target_Ref, fill=Strain)) +
geom_boxplot(alpha=0.6, outlier.shape = NA)+
geom_jitter(position=position_jitter(0.5), aes(color=Strain), alpha=0.6) +
theme_bw(base_size = 14)+
scale_fill_manual(values = c("blue", "grey", "red"))+
scale_color_manual(values = c("blue", "grey", "red"))+
labs(x="\nStrain", y="Ratio Target Ref (qPCR)", fill="Strain")+
guides(colour=FALSE)+
#stat_summary(fun = "median", geom = "point", shape = 23, size = 3, fill = "white")+
stat_compare_means(comparisons = my_comparisons,
method="wilcox.test",
label="p.signif",
vjust=1.5,
label.y = c(10, 11, 12))+
stat_compare_means(comparisons = my_comparisons,
method="wilcox.test",
label="p.format",
label.y = c(10, 11, 12))+
geom_text(data = means, aes(label = qPCR_Ratio_Target_Ref, y = qPCR_Ratio_Target_Ref + 0.4))
p8x <- sam3$Wolbachia %>% log10
y <- sam3$qPCR_Ratio_Target_Ref %>% log10
plot(x, y, pch = 19, col = "lightblue")
abline(lm(y ~ x), col = "red", lwd = 3)
text(paste("Correlation:", round(cor(x, y), 2)), x = -0.6, y = 0.2)# panels
p_group <- plot_grid(p2,
p3,
p4,
p5,
nrow=2,
ncol=2,
align="hv",
labels=c("A","B","C", "D")
)
p_group#ggarrange(p2, p3, p4, p5,nrow=2)
# save
setwd(path_plot)
tiff("1Ec_NMDS_main.tiff", units="in", width=8, height=5, res=300)
p6
dev.off()## quartz_off_screen
## 2
## quartz_off_screen
## 2
## quartz_off_screen
## 2
## quartz_off_screen
## 2
## quartz_off_screen
## 2
## quartz_off_screen
## 2
qPCR data
# all info
p1 <- plot_ordination(ps_percent_whole, bray.culex_whole, color="qPCR_Infection", shape=NA) +
labs(x="NMDS1", y = "NMDS2")+
stat_ellipse(geom = "polygon", level=0.95, alpha = 0.25, aes(fill = Species, color=NA))+
scale_fill_manual(values=c("orange", "yellow","green"), labels = c(expression(italic("Aedes aegypti")), expression(italic("Culex pipiens")), expression(italic("Culex quinquefasciatus"))))+
scale_color_manual(values=c("black", "blue", "red"), breaks=c("No infection", "Low infection", "High infection"))+
theme_gray()+
labs(caption=paste("stress:", bray.culex_whole$stress %>% round(2), "\n", "ellipse level: 0.95"))+
expand_limits(x=c(-2.5, 1.5), y=c(-2, 2.5))+
guides(fill = guide_legend(order = 1),
colour = guide_legend(order = 2),
shape = guide_legend(order = 3))+
labs(fill="Species", col="Wolbachia Group")+
geom_text(mapping = aes(label = Sample), size = 1.5, vjust = 3)+
theme(legend.text.align = 0)## Warning in plot_ordination(ps_percent_whole, bray.culex_whole, color =
## "qPCR_Infection", : Shape variable was not found in the available data you
## provided.No shape mapped.
p1$layers[[1]] <- NULL
p1[["guides"]][["shape"]][["title"]] <- "Strain"
p1[["guides"]][["colour"]][["title"]] <- expression(paste(italic("Wolbachia"), "Infection"))
p1+
geom_point(aes(shape = Strain), size = 2)## Warning: Removed 34 rows containing missing values (geom_point).
Statistical analysis
Effect of field
require(vegan)
setwd(path_xlsx)
# Select pipiens whole from
ps_pipiens <- ps.filter %>% subset_samples(Species=="Culex pipiens")
ps_pipiens <- check_ps(ps_pipiens)
ps_pipiens_whole <- ps_pipiens %>% subset_samples(Organ=="Whole")
ps_pipiens_whole_bosc_lavar <- subset_samples(ps_pipiens_whole, Strain!="Field - Camping Europe")
ps_pipiens_whole_camping_lavar <- subset_samples(ps_pipiens_whole, Strain!="Field - Bosc")
ps_pipiens_whole_bosc_camping <- subset_samples(ps_pipiens_whole, Strain!="Laboratory - Lavar")
# Adonis method
## Bosc / Lavar
adonis_field_bosc_lavar <- adonis(vegdist(t(otu_table(ps_pipiens_whole_bosc_lavar)), method = "bray") ~Field,
data=as(sample_data(ps_pipiens_whole_bosc_lavar), "data.frame"), permutation = 9999)
res_adonis1 <- adonis_field_bosc_lavar$aov.tab
res_adonis1 <- fill_significance(res_adonis1,"Pr(>F)")
## Camping / Lavar
adonis_field_camping_lavar <- adonis(vegdist(t(otu_table(ps_pipiens_whole_camping_lavar)), method = "bray") ~Field,
data=as(sample_data(ps_pipiens_whole_camping_lavar), "data.frame"), permutation = 9999)
res_adonis2 <- adonis_field_camping_lavar$aov.tab
res_adonis2 <- fill_significance(res_adonis2,"Pr(>F)")
## Bosc / Camping
adonis_field_bosc_camping <- adonis(vegdist(t(otu_table(ps_pipiens_whole_bosc_camping)), method = "bray") ~Strain,
data=as(sample_data(ps_pipiens_whole_bosc_camping), "data.frame"), permutation = 9999)
res_adonis3 <- adonis_field_bosc_camping$aov.tab
res_adonis3 <- fill_significance(res_adonis3,"Pr(>F)")
# Excel
wb <- loadWorkbook("Supplementary_Table_2.xlsx")
addWorksheet(wb, "NMDS")
writeData(wb, "NMDS", "Effect of lab (and Strain) in the Culex pipiens samples", startCol = 1, startRow = 25)
writeData(wb, "NMDS", "ADONIS (Analysis of variance using distance matrices):", startCol = 1, startRow = 27)
writeData(wb, "NMDS", "Lavar / Bosc:", startCol = 1, startRow = 29)
writeData(wb, "NMDS", res_adonis1, startCol = 1, startRow = 30, rowNames=TRUE)
writeData(wb, "NMDS", "Lavar / Camping Europe:", startCol = 1, startRow = 35)
writeData(wb, "NMDS", res_adonis2, startCol = 1, startRow = 36, rowNames=TRUE)
writeData(wb, "NMDS", "Bosc / Camping Europe:", startCol = 1, startRow = 41)
writeData(wb, "NMDS", res_adonis3, startCol = 1, startRow = 42, rowNames=TRUE)
saveWorkbook(wb, "Supplementary_Table_2.xlsx", overwrite = TRUE)Effect of tetracycline
setwd(path_xlsx)
# Select quinque whole
ps_quinque <- ps.filter %>% subset_samples(Species=="Culex quinquefasciatus")
ps_quinque <- check_ps(ps_quinque)
ps_quinque_whole <- ps_quinque %>% subset_samples(Organ=="Whole")
# Adonis method
adonis_antibiotics <- adonis(vegdist(t(otu_table(ps_quinque_whole)), method = "bray") ~Field,
data=as(sample_data(ps_quinque_whole), "data.frame"), permutation = 9999)
adonis_antibiotics$aov.tab ## Permutation: free
## Number of permutations: 9999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Field 1 1.4603 1.46029 5.191 0.25709 4e-04 ***
## Residuals 15 4.2197 0.28131 0.74291
## Total 16 5.6800 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
res_adonis5 <- adonis_antibiotics$aov.tab %>% as.data.frame()
res_adonis5 <- fill_significance(res_adonis3,"Pr(>F)")
writeData(wb, "NMDS", "Effect of tetracycline in the Culex quinquefasciatus samples", startCol = 1, startRow = 13)
writeData(wb, "NMDS", "ADONIS (Analysis of variance using distance matrices):", startCol = 1, startRow = 15)
writeData(wb, "NMDS", res_adonis3, startCol = 1, startRow = 17, rowNames=TRUE)
saveWorkbook(wb, "Supplementary_Table_2.xlsx", overwrite = TRUE)Effect of species
setwd(path_xlsx)
ps_whole_field <- ps.filter %>% subset_samples(Field=="Field")
adonis_whole_field <- adonis(vegdist(t(otu_table(ps_whole_field)), method = "bray") ~Species,
data=as(sample_data(ps_whole_field), "data.frame"), permutation = 9999)
res_adonis6 <- adonis_whole_field$aov.tab
res_adonis6 <- fill_significance(res_adonis6,"Pr(>F)")
writeData(wb, "NMDS", "Effect of species in the field samples", startCol = 1, startRow = 1)
writeData(wb, "NMDS", "ADONIS (Analysis of variance using distance matrices):", startCol = 1, startRow = 3)
writeData(wb, "NMDS", res_adonis6, startCol = 1, startRow = 5, rowNames=TRUE)
saveWorkbook(wb, "Supplementary_Table_2.xlsx", overwrite = TRUE)Effect of Wolbachia infection within Culex pipiens
setwd(path_xlsx)
# HCA Wolbachia infection * Strain
adonis_whole_infection <- adonis(vegdist(t(otu_table(ps_pipiens_whole)), method = "bray") ~Infection*Strain,
data=as(sample_data(ps_pipiens_whole), "data.frame"), permutation = 9999)
res_adonis7 <- adonis_whole_infection$aov.tab
res_adonis7 <- fill_significance(res_adonis7,"Pr(>F)")
writeData(wb, "NMDS", "Effect of Wolbachia infection (HCA groups) and Strain in Culex pipiens", startCol = 1, startRow = 50)
writeData(wb, "NMDS", "ADONIS (Analysis of variance using distance matrices):", startCol = 1, startRow = 52)
writeData(wb, "NMDS", res_adonis7, startCol = 1, startRow = 54, rowNames=TRUE)
saveWorkbook(wb, "Supplementary_Table_2.xlsx", overwrite = TRUE)Effect of Wolbachia infection, qPCR values and Strain within Culex pipiens with qPCR values
setwd(path_xlsx)
# HCA Wolbachia infection * Strain in Culex pipiens with a qPCR value
ps_pipiens_whole_qpcr <- ps_pipiens_whole %>% subset_samples(qPCR_Ratio_Target_Ref >0)
adonis_whole_infection2 <- adonis(vegdist(t(otu_table(ps_pipiens_whole_qpcr)), method = "bray") ~qPCR_Ratio_Target_Ref*Infection*Strain,
data=as(sample_data(ps_pipiens_whole_qpcr), "data.frame"), permutation = 9999)
res_adonis8 <- adonis_whole_infection2$aov.tab
res_adonis8 <- fill_significance(res_adonis8,"Pr(>F)")
writeData(wb, "NMDS", "Effect of Wolbachia infection (HCA groups), qPCR values and Strain in Culex pipiens with qPCR values", startCol = 1, startRow = 63)
writeData(wb, "NMDS", "ADONIS (Analysis of variance using distance matrices):", startCol = 1, startRow = 65)
writeData(wb, "NMDS", res_adonis8, startCol = 1, startRow = 67, rowNames=TRUE)
saveWorkbook(wb, "Supplementary_Table_2.xlsx", overwrite = TRUE)